Importing MNIST data

(* MNIST-images and MNIST-labels files in the notebook directory are required! *)

SetDirectory@NotebookDirectory[];

(* parsing the image file *)
rawimgs=ImportByteArray[ReadByteArray["MNIST-images"]];
read32[data_,off_]:=data[[#+1;;#+4]].Reverse[Power[2,8(Range[4]-1)]]&[off];
{mn,npics,nrows,ncols}=read32[rawimgs,#]&/@Range[0,12,4];
in=Partition[rawimgs[[17;;]]/255., nrows ncols]; (* converting the training images into input vectors *)
imgs=Image[Partition[#,ncols]]&/@in;

(* parsing the label file *)
rawlabs=ImportByteArray[ReadByteArray["MNIST-labels"]];
{mnl,nlabs}=read32[rawlabs,#]&/@Range[0,4,4];
labels=rawlabs[[9;;]];

(* previewing imported images and labels *)
imgs[[;;30]]//GraphicsRow
labels[[;;30]]

Neural Network parameters

(* Activation function (Sigmoid) *)
σ[x_]:=1/(1.+Exp[-x]);
(* Activation function derivative *)
dσ[x_]:=σ[x]/(1+Exp[x]);
(* Functions for each layer *)
f=σ;df=dσ; 
g=σ;dg=dσ;

(* Number of neurons for each layer *)
Ni=nrows ncols; (* input (image vector) *)
N1=28; (* hidden layer *)
N0=10;  (* output (one-hot encoding of a digit) *)

(* Initializing weights and biases *)
u=Array[RandomReal[{-1,1}]&,N0];
v=Array[RandomReal[{-1,1}]&,N1];
U=Array[RandomReal[{-1,1}]&,{N0,N1}];
V=Array[RandomReal[{-1,1}]&,{N1,Ni}];

Neural Network Training

(* Metaparameters *)
batchsize=100; 
α=5.; (* gradient descent speed *)

batches=Partition[{in,labels}\[Transpose],batchsize]; (* image(vector)/label pairs *)
dt=1./Length[batches]; 

(* calculating error gradient *)
BackProp[{u_,v_,U_,V_}]:=Function[{i,λ},Module[{dedu,dedv,dedU,dedV,dedb,l,ol,y,b,z,o},

(* Forward Pass *)
y=V.i+v;
b=g[y];
z=U.b+u;
o=f[z];

l=IdentityMatrix[10][[λ+1]]; (* converting label into one-hot encoding *)

(* Back Propagation *)
ol=(o-l);
dedu=ol*df[z]; (* ε=(o-l).(o-l); dedu=dedz; dodz=dfdz=df *)
dedU=dedu ⊗ b;
dedb=U\[Transpose].dedu;
dedv =dedb*dg[y]; (* dedv=dedy; dbdy=dgdy=dg *)
dedV=dedv ⊗ i;

{{dedu,dedv,dedU,dedV},ol.ol} (* output: {((∂m)/(δε)), 2ε} *)
]];

(* one step of stochastic gradient descent *)
GradStep[uvUV_,batch_]:=Module[{grade,res},

grade=Total[BackProp[uvUV]@@@batch]/batchsize; (* error gradient averaged over a batch *)
epoch+=dt; (* training progress *)
AppendTo[elist,{epoch,grade[[2]]}]; (* Training log *)
uvUV-α grade[[1]] (* updating weights and biases *)

];

(* Initializing training *)
epoch=0;
maxEpoch=3; (* training duration (number of full cycles over the training data) *)
elist={{0,1}}; (* initializing log *)

(* Training *)
Print[Dynamic[epoch]]
Dynamic[ListLinePlot[elist,PlotStyle->Black,PlotRange->{0,1}],SaveDefinitions->True] (* Training log graph *)
model=Nest[Fold[GradStep,#,batches]&,{u,v,U,V},maxEpoch]; (* Main loop *)

Testing training results

(* functions from the previous lesson *)
ForwPass[i_]:=Module[{u,v,U,V,y,b,z,o},
{u,v,U,V}=model;
y=V.i+v;
b=g[y];
z=U.b+u;
f[z]
];

OCR[img_]:=Module[{x,out},
x=Flatten@ImageData@ColorConvert[img,"Grayscale"];
out=ForwPass[x];
First@Ordering[out,-1]-1 
];

(* functions for graphical representation of recognition statistics *)
OCRcell=Graphics[{
White,Opacity[#[[2]]],Text[Style[#[[1]],FontSize->50,FontFamily->"Arial"],{0,.15}],
Opacity[.5],Text[Style[ToString@Round[100 #[[2]]]<>"%",FontSize->14,FontFamily->"Arial"],{0,-.65}]
},PlotRange->{{-1,1},1.1{-1,1}},Background->Black]&;

OCRstat[img_]:=GraphicsRow[Prepend[OCRcell/@SortBy[{Range[0,9],((#)/(Total[#]))&@ForwPass[img//ImageData//Flatten] }\[Transpose],-#[[2]]&][[;;4]],img],Spacings->0,ImageSize->300,Background->Black];
OCRstat[imgs__]:=OCRstat/@List[imgs];
OCRstat[imgs_List]:=OCRstat/@imgs;

(* Recognition of the training images *)
OCRstat[{,,,,,,,,,}]//TableForm

(* Recognition of new images *)
OCRstat[{,,,,,,,,,}]//TableForm